function [pVentral,pDorsal] = initialPosition(totAngle,lh,lv,N)
% function p = initialPosition(totAngle,lh,N)
% totAngle in radians
%vector of lh should have N-1 elements
%vector of lv should have N elements

L = sum(lh);
R = L / totAngle; 

deltaTheta(2:N) = 2*asin(lh(1:N-1)/(2*R));
deltaTheta(1) = 0;
theta = zeros(1,N);

theta(1:N) = cumsum(deltaTheta);

x = R*cos(theta-(pi/2));
y = R*sin(theta-(pi/2))+R;

xMinus(1:N) = x(1:N) + 0.5*lv(1:N) .* cos ( theta-(pi/2)); 
yMinus(1:N) = y(1:N) + 0.5*lv(1:N).* sin ( theta-(pi/2)) ;


xPlus(1:N) = x(1:N) - (lv(1:N) ./2).* cos ( theta-(pi/2)) ;
yPlus(1:N) = y(1:N) - (lv(1:N)./2).* sin ( theta-(pi/2)) ;


pDorsal=[xPlus;yPlus]';
pVentral = [xMinus;yMinus]';


%Plotting the initial position for testing
%plot(x,y,'-or'); hold on
%plot(xPlus,yPlus,'-og'); 
%plot(xMinus,yMinus,'-ob'); hold off;
%hold off;
 %axis square